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I.  PLASMA  THEORY  AND  SIMULATION 


A.  Lower  Hybrid  Drift  Instability;  2d  Simulations  and  Nonlocal  Theory. 

Yu-Jiuan  Chen  (Prof.  C.K.  Birdsall,  W.M.  Nevins,  LLL) 

The  lower  hybrid  drift  instability  was  studied  with  both  a  2d  electrostatic  simulation  code 
and  a  nonlocal  theory.  Results  will  be  given  in  an  ERL  report  now  in  preparation,  to  be  sub¬ 
mitted  to  the  Physics  of  Fluids. 

B.  Velocity  Space  Aliasing. 

Niels  Otani,  Dr.  B.I.  Cohen  (LLL),  and 
Dr.  M.J.  Gerver  (MIT)  (Prof.  C.K.  Birdsall) 

See  contribution  from  LLL,  following. 
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Velocity- Space  Aliasing  Instability  in  Magnetized  Quiet-Start  Simulations 
Bruce  I.  Cohen,  William  M.  Nevins,  and  Neil  Maron 
Lawrence  Livermore  National  Laboratory,  University  of  California 
Livermore,  California  94550 

In  a  continuing  effort  to  develop  a  magnetized  quiet -start  simulation 
technique,  we  have  studied  analytically  and  in  simulation  the  properties  of 

A 

distribution  functions  of  velocities  perpendicular  to  BQz  composed  of 
discrete  speed  classes  (rings)  and  angle  classes  (spokes).  A  previous 
progress  report1  offered  sufficiency  conditions  for  stability  of  equally 
weighted  rings  and  spokes  Maxwellian  distribution  functions.  Roughly 

o 

speaking,  if  the  number  of  velocity  rings  Ny^  satisfies  Ny^  >  (1/3) (u>p/u)c) 
and  the  number  of  velocity  spokes  N0  satisfies  NQ  >  8(kmaxvtf/a)c^  then 
stability  of  electrostatic  flute  modes  (k  •  B  *  0)  is  expected  (vth  =  thermal 
velocity,  u>p  =  plasma  frequency,  wc  =  cyclotron  frequency,  kmax  e  maximum 
wavenumber). 

The  stability  of  flute  modes  with  respect  to  the  number  of  velocity  rings 

2  2 

N  and  relative  density  u)p/u)£  was  numerically  determined  by  Kim  and  Otani 

2  1 

neglecting  the  effects  of  discrete  spokes.  Otani  and  Cohen  derived  a 
dispersion  relation  including  the  coupling  produced  by  discrete  spokes, 

47m  r  /  nd<J^/du 

‘  Jdvf o< “)  j *r  £  (n  -  r»c)  -  n»e 

♦e 

n,m 
m»0 


pr-mN, 


n3(JnJn-mNfl)/3y  "  mNeJn-mNe3jn/3w) 
V  _ I _ 


(u  -  r«c)  -  nwc  + 


mNe“c 
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(1) 


where  u  =mv2/2u.,  <Kx,t)  =  ^  4>r  exp(ikx  +  iw.t  -  iwt)  +  c.c., 

r 

and  the  Bessel  functions  have  argument  kvj/w  .  The  effective  frequency  of 

mode  is  ^  m  -  r^c .  Coupling  to  the  spoke-induced  cyclotron  aliases  will 

be  negligible  if  the  product  is  small. 

The  coupling  of  angular  spokes  deserves  more  discussion.  For  large  index, 

the  Bessel  function  is  small  in  general  and  especially  small  for  argument 

kv,/u)c  less  than  the  absolute  value  of  the  index.  For  argument  much  greater 

-1/2 

than  the  index,  the  magnitude  of  the  Bessel  function  decreases  as  (kv^/u^)  ' 

and  the  magnitude  of  the  Bessel  function  is  maximized  for  kvj_/w  comparable 
to  the  absolute  value  of  the  index.  Possible  couplings  to  cyclotron  harmonics 
separated  by  +  N0u>c(m  ■  +  1)  are  the  most  likely  to  occur,  if  any  alias 
coupling  is  going  to  occur  at  all.  The  arguments  of  the  Bessel  functions  are 
naturally  limited  by  the  largest  wavenumber  and  velocity 


kv./w.  <  4k  v.. /u). 
1  c  max  th  c 


(2) 


The  product  JnJn+N  will  be  small  for  all  values  of  kVj^/a^if 
—  0 

4kmaxvth/uc  <min  [lnl»  In-Ne0  1  V2  *  (3) 

> 

If  Eq.  (3)  is  not  satisfied,  cyclotron  alias  coupling  can  occur  and 

instability  may  result.  The  first  modes  to  couple  as  we  increase  the  left 

side  of  Eq.  (3)  above  the  stability  limit  satisfy  |n|  *  iN^+n |  *  N0/2 ;  hence, 

w  as  +  N jji/2.  A  simulation  example  Is  pictured  in  Fig.  1.  In  this 
r  —  o 

simulation,  the  parameters  were  o^/u2  ■  200,  NQ  -  32,  and  *  299.  The 
distribution  function  was  replicated  eight  times  across  the  length  of  the 


system  to  quiet  the  first  few  spatial  modes.  The  numerical  solutions  of  Kim 


6 


2 

and  Otani  suggest  that  the  distribution  of  rings  should  have  been  stable. 

For  the  mode  pictured  in  Fig.  1,  =  6.25  and  velocities  in  the  tail 

of  the  distribution  2.5  v^  ±  v^  £  3  v^  give  arguments  falling  near  the 
maxima  of  J^,  for  14  |«.|  <,  18.  Use  of  the  ZED  post-processor  revealed  the 

signature  of  the  cyclotron  alias  coupling  as  lower  hybrid  normal  modes  near  * 

(?«“*  +  14“  coupled  to  modes  at  Re  “  =  +  (14  -  32)“_  =  +  18 

i  C  *  C  C  0 

these  modes  all  grew  with  rate  Im“r  «  0.15  “c  and  the  ratios  of  the 

coupling  modes  were  equal  as  expected,  |<{>(14  “  ,k)p/l«l>(-18  “c,k)r  = 
l*(-14  “c,k)l2/|  d>(18  “c,k) 1 2. 

We  doubled  NQ  to  NQ  =64,  so  that  4kmaxvth/aJc  =  25,  NQ/2  =  32,  and 

Eq.  (3)  was  satisfied.  Figure  2  shows  that  modes  near  the  lower  hybrid 

frequencies  Re  w  =  +  14  u  were  excited  and  may  have  been  very  weakly 

unstable.  However,  if  there  was  instability,  saturation  occurred  very  quickly 

2  -11 

and  at  a  very  small  amplitude  E  /8jrnT  =  10  ,  which  was  not  much  above 

the  initial  amplitude.  In  fact,  the  fluctuation  level  in  the  normal  modes  was 
so  small  that  the  discrete  particle  ballistic  noise  was  visible  in  the  power 
spectrum  (Fig.  2b).  For  our  purposes,  this  simulation  corresponded  to  a 
stable  quiet  start.  We  note  in  closing  that  Homer  Meier  explored  some  of 
these  issues  ten  years  ago. 


♦This  work  was  performed  under  the  auspices  of  the  U.S.  Department  of  Energy 
at  the  Lawrence  Livermore  National  Laboratory  under  contract  W-7405-ENG-48. 

1.  N.  F.  Otani  and  B.  I.  Cohen,  "An  Electrostatic  Dispersion  Relation  for  e 
Uniform  Magnetized  Plasma  Having  a  Discrete-Particle  Velocity 
Distribution,"  Fourth  Quarter  Progress  Report  on  Plasma  Theory  and 
Simulation,  October  1  to  December  31,  1980,  Electronics  Research 
Laboratory,  U.C.  Berkeley. 


2.  J.  S.  Kim  and  N.  Otani,  "Magnetized  Multi-Ring  Instabilities,"  Second 


Quarter  Progress  Report  on  Plasma  Theory  and  Simulation,  April  1  to  June 
30,  1980,  Electronics  Research  Laboratory,  U.C.  Berkeley. 

3.  H.  Meier,  unpublished. 


Figure  Captions 

Figure  1.  Simulation  of  a  Maxwellian  plasma  with  a  separable  spokes  and  rings 
perpendicular  velocity  distribution  function  for  electrostatic  flute  modes 
(k  *  B  =  0).  The  parameters  were  =  200,  kvth/u>c  =  1.25,2.5,3.75,5, 

and  6.25,  Nv  =  299  and  N0  -  32.  (a)  Perpendicular  velocity  distribution 

function.  Apparent  structure  at  low  velocities  is  due  to  the  particular 
sampling  of  3000  particles  out  of  the  76544  total,  (b)  Power  spectrum 

O 

|4>(oj,k)|  as  a  function  of  the  frequency  ui  and  (c)  The  corresponding  field 
energy  E£/8it  integrated  over  u  as  a  function  of  time  <i>ct  for  kvth/o)c  =  6.25. 

Figure  2.  Effectively  stable  spokes  and  rings  simulation  with  the  same  plasma 
parameters  as  in  Fig.  1  and  with  N  =  149,  NQ  =  64  and  76288  total 

Vj^  y 

particles,  (a)  Sampling  of  perpendicular  velocity  distribution  function,  (b) 

Power  spectrum  | <t> (co,k ) |  as  a  function  of  frequency  w  and  (c)  E£/8ir  as 

a  function  of  time  for  kv«.  =  6.25. 

c  th  c 
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C.  Hybrid  Simulations  of  Theta  Pinch  Rotational  Instabilities. 

Douglas  Harned  (Prof.  C.  K.  Birdsall) 

Rotating  rigid  rotor  theta  pinch  equilibria  have  been  initialized  in  our  two-dimensional 
hybrid  code,  AQUARIUS.  The  code  now  performs  a  three  region  solution  of  the  field  equa¬ 
tions  which  properly  treat  the  plasma-vacuum  interface.  Stationary  equilibria  have  been  found 
to  be  stable  to  all  rotational  instabilities,  while  equilibria  with  co-rotating  ions  and  electrons 
have  been  found  to  develop  an  unstable  m— 2  mode.  Growth  rates  have  been  measured  and 
appear  to  be  comparable  to  those  predicted  by  the  Vlasov- fluid  theory  of  Seyler1  and  the  finite 
Larmor  radius  theory  of  Freidberg  and  Pearlstein2. 

As  the  m—2  mode  grows  to  large  amplitude,  nonlinear  effects  reduce  the  growth  rate. 
This  appears  to  be  primarily  due  to  stabilization  by  the  conducting  wall.  However,  the  instabil¬ 
ity  continues  to  grow  at  a  reduced  rate  and  leads  to  the  eventual  disruption  of  the  plasma. 

References 

1.  C.E.  Seyler,  "Vlasov-Fluid  Stabilily  of  a  Rigidly  Rotating  Theta  Pinch,"  Phys.  Fluids  22,  2324,  (1979). 

2.  J.P.  Freidberg,  L.D.  Pearlstein.  "Rotational  Instabilities  in  a  Theta  Pinch,"  Phys.  Fluids  21,  1207,  (1978). 


II.  CODE  DEVELOPMENT  AND  MAINTENANCE 
A.  Implicit  Particle  Simulation  using  Moments  with  Orbit  Averaging. 
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V.  Thomas  (Dr.  B.  1.  Cohen  LLL  and  Prof.  C.  K.  Birdsall) 


In  our  previous  QPR  we  reported  our  initial  failure  to  obtain  a  stable  orbit  averaged  impli¬ 
cit  moment  code.  We  have  since  implemented  many  other  algorithms  in  an  attempt  to  create  a 
stable  code.  One  stable  algorithm  was  found  so  far.  We  have  developed  model  equations  for 
our  algorithms  which  correctly  predict  their  stability  for  the  case  of  linear  cold  plasma  oscilla¬ 
tions.  Our  stable  orbit  averaged  code  was  used  successfully  to  simulate  a  single  mode  ion 
acoustic  wave.  However,  although  we  have  constructed  a  stable  orbit  averaged  implicit  moment 
code,  we  apparently  obtain  no  advantages  over  the  non-orbit  averaged  implicit  moment  code  as 
used  by  Mason  for  electrostatic  wave  propagation  parallel  to  the  magnetic  field.  The  orbit  aver¬ 
aged  code  is  smaller  but  the  number  of  operations  cannot  be  reduced.  We  note  in  passing  that 
for  wave  propagation  parallel  to  the  magnetic  field  a  simple  lag  average  of  the  pressure  and 
current  from  previous  time  steps  may  prove  effective  in  reducing  the  number  of  particles 
required  for  an  implicit  moment  code.  We  have  tried  this  in  an  implicit  moment  code  and  it  is 
stable  (except  for  a  lag  average  of  the  density,  which  is  apparently  always  unstable  ).  Further 
study  is  needed  to  test  for  savings. 


The  stable  orbit  averaged  algorithm  employs  particle  ions  and  electrons  to  calculate  the 
current  density  and  the  kinetic  stress  tensor  on  the  simulation  grid.  The  charge  density  is  a 
fluid  quantity  and  it  is  never  accumulated  from  the  particles  in  a  direct  fashion. 


At  a  given  time  step  A,  the  current  density  j  is  known  at  time  step  A- 1/2,  the  kinetic 
stress  tensor  is  known  at  time  A- 1/2,  and  the  fluid  density  is  known  at  time  step  A.  The 
current  density  and  the  kinetic  stress  tensor  are  orbit  averaged  quantities  and  they  are  given  by 


iN-\/2  - 


1 


mjv+i>AI 


pN- 1/2  = 

*i  — 


(i  +  AI) 

A/s  H 

ar 

a  i. 


/-( JV+l) 


1 


(1  +  AI)  ,-*A L 

^ts  A', 


Pi 


(i) 


(2) 


Here  s  is  the  species  label,  A  T  is  the  field  solve  time  interval,  and  A  t5  is  the  micro  time  step 
for  species  s. 

The  fluid  density  at  time  A+l  is  obtained  implicitly  in  terms  of  EN+X  through  the  use  of 


\+ 1 
S 


esn 
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#-y*+,/JAr 

dx 


jrm 


JLp/V-l/2 
dx  s 


A  T 


(3) 

(4) 


Here  E'  represents  a  linear  combination  of  EN+\  EN  ,  and  EN  1  ,  Js  represents  the  fluid 
current  density,  and represents  the  particle  current  density. 

Next  Poisson’s  equation  is  used  to  solve  for  the  implicit  field  £,'V+I  as  in 


EN*l’-f'Zesn?+'dx+  C  (5) 

0  j 

Here  C  is  a  constant  to  be  determined  from  the  boundary  conditions.  The  particles  are  then 
advanced  and  the  pressure  term  is  recalculated.  Equations  (3)  and  (4)  are  used  again  to  deter¬ 
mine  a  new  n?+x  and  EN+K  The  particles  are  then  re-advanced  and  the  algorithm  proceeds  to 


the  next  time  step. 

The  model  equations  for  the  stable  orbit  averaging  simulation  method  are  given  in  the 
following  four  equations.  The  equations  are  valid  for  cold  linear  plasma  oscillations  in  the  limit 
of  1 tAx  —  0. 


i/c£M+ 1  -  ebnN+l 
e(,&nN+l  -  8nN)  -  -ikJN+y2AT 


;W./-W+f!54r£(.  AVN~'  +  v")  +  e2— A7£* 
m2  m 

v*+1  -  vN  +  —A  TE* 
m 

£*  -  a£*+l  +  0 En  +  (!-«-  p)EN~l 


(6) 

(7) 

(8) 

(9) 

(10) 


The  dispersion  relation  we  obtain  is 

-( z  ~  l)2  -  at  jA  r2|z  +  -  y  az  +  p+  #  (11) 

Here  z  -  e  ~  '"4  r  and  a  and  p  are  parameters. 

Two  test  runs  were  made  for  the  stable  algorithm.  The  amplification  factors  obtained 
from  the  dispersion  relation  agreed  with  those  obtained  by  simulation  to  three  places.  For 
aip A  T  -  0.2  the  simulation  value  for  the  amplification  factor  was  0.990  whereas  for  the  disper¬ 
sion  relation  we  expect  the  value  0.990.  For  <upA  7—4,  the  simulation  value  for  the 
amplification  factor  was  0.727  versus  the  value  predicted  by  the  dispersion  relation  of  0.728. 
Comparison  of  the  model  equations  with  3  unstable  schemes  gave  results  almost  as  good  as 
thesf,.  All  simulations  were  performed  using  infinitely  massive  ions  and  cold  electrons  with 
excitations  of  approximately  10~4of  a  particle  separation  distance. 

One  concern  for  the  stable  algorithm  is  the  possible  divergence  of  the  fluid  density  with 
the  particle  density,  since  the  fluid  density  is  never  accumulated  from  the  particles.  To  test  this 
over  a  long  time  scale  we  performed  a  single  mode  ion  acoustic  simulation  to  the  time  24000 
(Op  The  fluid  charge  density  and  the  particle  charge  density  accumulated  form  the  grid  are 
shown  in  Figure  1  and  in  Figure  2  respectively.  The  results  appear  to  be  satisfactory.  The  plot 
of  the  ion  kinetic  energy  vs.  time  in  Figure  3  exhibits  the  correct  frequency.  The  theoretical 
damping  rate  is  not  easily  seen  from  the  plot.  The  field  energy  for  this  run  is  shown  in  Figure 
4. 

For  the  ion  acoustic  runs  that  used  E  equal  to  £'v+l,  the  algorithm  was  unstable  for 
*max  V»A  T  ^  0.1  in  contrast  to  the  method  employed  by  Mason  where  there  is  no  such  stability 
limit  when  a  fully  implicit  field  corresponding  to  our  EN+l  is  used.  With  further  modification 
of  the  orbit  averaging  algorithm  it  might  be  possible  to  relax  the  v* A  T  constraint. 

We  have  simulated  a  single  mode  ion  acoustic  wave  using  Mason’s  method  to  compare 
with  our  orbit  averaged  simulation.  Figure  5  shows  the  ion  kinetic  energy  for  the  simulation 
and  Figure  6  shows  the  field  energy. 

We  chose  the  number  of  particles  for  this  simulation  so  that  the  orbit  averaged  simulation 
and  this  simulation  would  have  the  same  number  of  particle  pushes  (i.e.  so  that  both  simula¬ 
tions  would  require  a  similar  number  of  computer  operations,  ignoring  for  the  moment  the  fact 
that  the  orbit  averaged  code  has  a  predictor  and  a  corrector  step  and  no  corrector  step  was  used 
for  our  Mason's  method  simulation).  From  Figures  3  to  6  it  is  clear  that  orbit  averaging  has 
not  helped  to  eliminate  noise  and  furthermore,  resonant  particle  effects  appear  to  be  better 
represented  by  our  Mason’s  method  simulation.  The  ion  kinetic  energy  from  our  Mason’s 
method  simulation  shows  damping  throughout  the  entire  simulation  from  t  -  0  to 
i  -  24000  Up 1  ,  although  the  damping  is  slower  than  the  Landau  rate.  This  may  be  due  to 
finite  grid  effects.  In  contrast  there  is  no  damping  in  the  orbit  averaged  code.  In  fact  the  ion 
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Figure  1.  The  particle  charge  density  averaged  over  a  macrostep  at  time  /  -  2O8OO0), For  this  simulation  we  have 
used;  mj m,  —  25,  T  —  4,  Af,  —  A/,  -  0.2,  v„  —  0.025,  v#  —  0.0,  NG  —  32,  Ax  -  1,  /Vp  -  512,  £*—  £A,+I,  and 
the  ions  were  given  an  initial  sinusoidal  velocity  perturbation  in  the  simulation  mode  (mode  I)  with  amplitude 
V|  -  -0.0005.  The  fluid  density  was  initialized  in  a  self  consistent  manner. 


0  8  16  24  32 
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Figure  2.  The  fluid  charge  density  at  time  /  -  20800aij'. 

kinetic  energy  reaches  a  maximum  of  approximately  50%  larger  than  the  initial  excitation  at  the 
end  of  the  orbit  averaged  simulation.  Also  the  noise  in  the  held  energy  is  greater  for  the  orbit 
averaged  simulation  than  for  the  Masons’  method  simulation.  This  is  due  to  the  fact  that  the 
particles  are  not  able  to  hi!  in  much  of  phase  space  during  a  macrostep  because  of  the  very  res¬ 
trictive  condition  on  T. 

Orbit  averaging  cannot  improve  particle  statistics  over  the  implicit  moment  particle  simu¬ 
lation  method  in  the  unmagnetized  electrostatic  limit  because  the  moment  implicit  method  uses 
a  held  solve  time  step  that  is  the  same  as  the  desired  particle  time  step;  i.e.,  the  largest  A  T  that 
still  accurately  resolves  the  particle  orbit.  The  source  term  for  the  orbit  average  electric  held 
contains  (A  T/£kt)NP\  input  variables,  where  NP\  is  the  number  of  particles.  The  source  term 
for  the  implicit  moment  simulation  method  contains  NP2  input  variables.  In  order  to  achieve 
the  same  statistical  behavior  in  the  two  sources  we  might  expect  to  require  the  same  number  of 
particle  pushes,  T/ A bfP  1  —  NP2  ,  in  which  case  averaging  vjjild  not  give  savings,  even  for 
the  same  A  T.  The  situation  with  this  orbit  averaged  code  is,  in  fact,  worse  than  this  because 
one  may  use  a  larger  A  T  with  the  implicit  moment  method  and  still  retain  stability  (although 
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Figure  3.  The  ion  kinetic  "sloshing"  energy  from  /  -  0  to  r  —  1 2000^ ~ 1 


Figure  4.  The  electrostatic  energy  from  i  -  0  to  r  —  12000u^'. 

we  still  require  kmtxv„ii  T  «  1  to  retain  all  of  the  physics  ). 

These  results  are  in  contrast  to  the  results  obtained  when  orbit  averaging  was  used  in  a 
magneto-inductive  code.*  In  that  instance  the  number  of  particles  required  for  the  simulation 
was  dramatically  reduced  due  to  the  fact  that  A  T  was  made  very  large.  Also,  since  the  particles 
were  moved  with  ft  fA  t  «  1  the  cyclotron  motion  was  treated  accurately. 

For  cases  where  one  needs  to  accurately  treat  wave  phenomena  parallel  to  a  magnetic  field 
or  in  its  absence,  one  needs  k^  v*A  T  «  1.  If.  this  is  the  most  severe  constraint  on  if  for 
accurate  particle  motion,  it  is  impossible  for  orbit  averaging  to  give  any  great  advantage  if  one 
can  move  the  particles  in  time  steps  of  A  T  with  the  use  of  another  method.  For  cases  where 
flrA7’  »  1,  orbit  averaging  gives  more  accurate  gyromotion. 

It  is  important 

to  note  that  these  remarks  are  true  when  the  wave  has  a  component  along  the  field  lines.  For 
the  other  case  orbit  averaging  can  offer  advantages.**  The  reader  is  referred  to  that  paper  for 


•Cohen  B.  I.,  Brengle  T.  A.,  Conley  D.  B.,  and  Freis  R.  P.,  *  An  Orbit  Averaged  Particle  Code,"  J.  Compul. 

Phys.  3B.  45,  Nov.  1980. 

••Freis  R.  P.,  Cohen  B.  I.,  Byers  J.  A.,  and  Thomas  V.  A.,  "Orbit  Averaged  and  Implicit  Particle  Codes", 

Annual  Sherwood  Controlled  Fusion  Theory  Conference  ,  paper  ICII  Univ.  of  Texas,  Austin,  Texas,  April  1981, 


Figure  5.  The  ion  kinetic  "  sloshing  ”  energy  from  t  -  0  to  i  -  8000m,, '.  Here,  Np  —  10240  and  the  other  parameters 
are  the  same  as  for  the  orbit  averaged  simulation. 


Figure  6.  The  electrostatic  energy  from  /  -  0  to  /  -  8000m 

that  case. 

B.  POLARES. 

Niels  F.  Otani  (C.  K.  Birdsall) 
No  special  progress  to  report. 

C.  PLASMA  SHEATH  SIMULATION  IN  Id. 

S.  Rousset  (Prof.  C.  K.  Birdsall) 


Modifications  to  the  ESI  code  have  been  written  so  as  to  allow  the  simulation  of  a  plasma 
sheath.  In  the  model  used,  we  have  a  plasma  reservoir  in  the  half-space  x<0  and  a  metal 
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target  at  x—L.  Depending  on  the  externally  applied  potential  between  x-0  and  x—L,  on  the 
external  magnetic  held  and  on  the  initial  velocity,  a  particle  can  or  cannot  reach  the  target.  The 
purpose  of  this  study  is  to  describe  the  velocity  distribution  on  the  target,  the  current  density, 
instabilities,  etc. 

In  order  to  have  a  periodic  numerical  model  to  describe  this  nonperiodic  physical  model, 
we  first  compute  the  true  charge  density  in  the  0<x<L  region,  and  we  add  an  inverted  charge 
density  in  the  L<x<2L  region  so  that  p(L+x)—p(L-x).  This  is  the  same  as  considering 
that  for  each  charge  in  the  true  sheath,  there  is  an  image  charge  (opposite  sign)  in  the  inverted 
sheath. 

Note  that  only  the  charge  density  and  the  fields  are  computed  in  the  inverted  region; 
there  is  actually  no  particle  outside  the  true  region:  the  equations  of  motion  are  integrated  in 
the  0<x<Z.  region  only,  and  whenever  the  mover  brings  a  particle  from  a  point  within  the 
sheath  to  a  point  x<0  or  x>  L,  the  particle  disappears. 

At  present  the  code  is  in  its  preliminary  testing  stages.  No  major  problem  has  yet  been 
encountered. 

References 
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D.  ICFT. 


H.  Stephen  Au-Yeung 

ICFT  is  a  preprocessor  of  CFT.  It  allows  the  user  to  specify  include  files  within  the 
source  program. 

An  include  file  contains  a  sequence  of  Fortran  statements.  For  those  who  are  familiar 
with  the  CHATR  compiler  on  the  CDC-7600,  an  include  file  has  similar  functions  to  a  cliche 
(see  the  CHATR  manual).  To  refer  to  an  include  file,  the  user  puts  the  include  statement  in 
one  of  the  two  following  formats: 

(1)  Within  a  comment  statement  - 
c  include  FileName 

(2)  Not  within  a  comment  statement  - 
include  FileName 

When  ICFT  reaches  an  include  statement,  it  opens  the  specified  file  and  replaces  the  statement 
by  whatever  is  in  the  file.  After  replacing  all  the  include  statements,  ICFT  then  invokes  CFT 
to  perform  the  compilation. 

ICFT  resides  in  the  FILEM  directory  .CRAY  of  user  number  1222.  To  obtain  ICFT  from 
the  CRAY-1  computer,  the  user  may  type: 

FILEM  READ  1222  CRAY  ICFT(ESC)END 
An  example  may  make  this  document  clearer: 

The  following  is  a  main  program  contained  in  the  file  •  SOURCE. 

c  An  example  to  show  how  to  use  the  include  file, 
c 

c  include  identify 


parameter  (  CPW-g,  BPW-64  ) 
parameter  (  NS  — 10- ) 


integer  Line (10) 
c 

include  initial 
c 

call  msgior  ( 'Example  of  include  file',  23  ) 

call  exit 

end 

The  file  IDENTIFY  contains 

c  AUTHOR:  H.  Stephen  Au-Yeung 

c  INSTITUTION.  University  of  California,  Berkeley 

and  the  file  INITIAL  contains 

call  dropfile  (  0  ) 
call  msglink  (  59,  1  ) 

The  input  on  the  execution  line  of  ICFT  is  the  same  as  that  of  CFT  plus  an  optional  ILIB 
entry.  A  typical  example  is: 

ICFT  I  -SOURCE.B  -  BIN  AR  Y,L  -  0,ILIB  -  (ICFT.INCLUDE) 

where  ICFT  and  INCLUDE  are  both  LIB  hies.  If  an  include  hie  is  not  found  under  the  active 
hie  index,  ICFT  opens  the  libraries  specified  in  the  ILIB  entry,  ICFT  and  INCLUDE  in  this 
case,  and  looks  for  the  hie.  If  there  needs  only  one  include  library,  the  parentheses  can  be 
omitted.  When  no  ILIB  entry  is  given,  only  the  active  hie  index  is  searched. 

Notes: 

(1)  Include  statements  may  also  appear  within  include  hies.  The  nesting  can  go  as  deep  as  10. 

(2)  Each  include  iibrary  may  contains  at  most  200  entries. 

E.  EMI  Code. 

No  special  progress  to  report. 

F.  ESI  Code. 

No  special  progress  to  report. 

G.  EZOHAR  Code. 


No  special  progress  to  report. 


18 

III.  SUMMARY  OF  REPORTS,  TALKS,  PUBLICATIONS,  VISITORS 

Journal  publication: 

Nevins,  W.M.,  Y.  Matsuda,  and  M.J.  Gerver,  ’Plasma  Simulations  Using  Inversion  Symmetry 
as  a  Boundary  Condition,'  J.  Comput.  Phys.  39,  1,  January  1981. 

Abstracts  of  3  talks  submitted  to  conferences  follow. 
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April  8-10,  1981 


'  Lower  Hybrid  Drift  Instability;  Theory  and  2d  Simulations 

Yu-Jiuan  Chen  and  C.  K.  Birds all 
Electronics  Research  Laboratory 
'University  of  California 
Berkeley,  CA  94720 

William  M.  Nevins 

Lawrence  Livermore  National  Laboratory 
University  of  California 
.  Livermore,  CA  94550 


The  lower  hybrid  drift  instability  is  being  studied  with  both  a  2-d 
electrostatic  simulation  code  and  a  nonlocal  theory.  A  slab  configuration  is  used 
with  a  density  gradient  in  x.  In  the  initial  equilibrium,  the  ion  pressure 
gradient  cancels  the  zeroth  order  ambipolar  electric  force.  Since  the  charac¬ 
teristic  frequency  of  the  lover  hybrid  drift  instability  is  much  greater  than 
the  ion  cyclotron  frequency,  the  ions  are  treated  as  unmagnetized  particles . 

Simulations  show  good  agreement  of  the  measured  growth  rate  and  frequency 
with  the  results  of  local  theory  during  the  early  stage  of  wave  growth.  At 
later  times  nonlocal  effects  become  important,  and  a  coherent  mode  structure 
deve lopes .  This  normal  mode  is  observed  to  propagate  up  the  density  gradient. 

The  mean  wave  vector  is  k-k  (-R  +?  ) ,  in  good  agreement  with  our  nonlocal 
theory .  ~  ^ 

At  zero  plasma  beta  and  zero  electron  temperature,  we  found  that  the  lower 
hybrid  drift  instability  is  stabilized  by  the  local  current  relaxation  occuring 
in  two  separate  regions.  The  neighborhood  of  the  point  where  the  electrons 
have  the  greatest  relative  drift  velocity  is  the  most  unstable  region  according 
to  local  theory.  In  this  region  the  current  relaxes  by  modification  of  the 
ion  velocity  distribution.  However,  in  our  nonlocal  theory,  the  wave  propagates 
into  the  plasma  to  a  second  region  where  the  electron  drift  velocity  equals 
the  phase  velocity  of  the  most  unstable  wave.  In  this  second  region  the  elec¬ 
trons  are  in  resonance  with  the  wave  and  local  current  relaxation  occurs  by 
modification  of  the  electron  density  profile. 


*  This  research  was  supported  in  part  by  the  Office  of  Naval  Research  under 
Contract  N00014-77-C-0578  (Berkeley) ,  and  in  part  by  the  Department  of 
Energy  under  Contract  No.  W-7405-ENG-48  (Livermore) . 

Submitted  for  poster  session,  Sherwood  Theory  Meeting,  Annual  Controlled 
Fusion  Theory  Conference,  Austin,  Texas,  April  8-10,  1981. 


Submitted  to  Sherwood  Theory  Meeting,  Austin,  'Ttexas,  April  8-10,  1981 

,  ROTATIONAL  INSTABILITIES  IN  FIELD-REVERSED  CONFIGURATIONS; 

RESULTS  OF  TWO-DIMENSIONAL  HYBRID  SIMULATIONS 

Douglas  S.  Harned 
Electronics  Research  Laboratory 
University  of  California,  Berkeley 


Rotational  instabilities  in  field-reversed  configurations 
are  being  studied  using  a  two-dimensional  hybrid  simulation  code 
The  code  treats  ions  as  particles  and  electrons  as  an  inertia¬ 
less  fluid.  Fields  are  described  by  the  Darwin  version  of  Max¬ 
well's  equations  (i.e.,  the  transverse  displacement  current  is 
neglected)  and  the  electron  momentum  equation,  coupled  by  the 
assumption  of  quasineutrality.  Ions  are  moved  by  the  equations 
of  motion.  Linear  weighting  is  used  to  determine  forces,  den¬ 
sities,  and  ion  currents.  The  code  has  been  applied  to  two  con¬ 
figurations!  (1)  a  tenuous  high-current  ion  layer  immersed  in 
a  dense  background  plasma,  and  (2)  a  field-reversed  theta-pinch. 

In  field-reversed  ion  layers  instabilities  may  arise  due 
to  a  coupling  of  layer  betatron  oscillations  with  compressional 
Alfven  waves  in  the  background  plasma.  The  hybrid  simulations, 
as  well  as  a  numerical  extension  of  previous  theoretical  work, 
have  demonstrated  that  instabilities  may  occur  in  thin  layers  at 
levels  of  field-reversal  below  those  previously  predicted.  This 
result  implies  that  field-reversed  ion  layers  must  be  very  thick 
to  be  8 table.  The  simulations  show  that  instability  growth 
causes  heating  within  the  layer,  which  results  in  the  layer  ex¬ 
panding,  reducing  the  betatron  frequency.  The  shift  in  betatron 
frequency  prevents  the  resonant  interaction  of  particles  with 
the  compressional  Alfven  waves,  ending  instability  growth.  Af¬ 
ter  all  instability  growth  has  stopped  the  layer  is  still  field- 
reversed,  but  much  thicker.  The  final  state  appears  to  be  rela¬ 
tively  stable  and  is  characterized  by  large  electron  currents 
and  many  non-axis- encircling  ions. 

The  field-reversed  theta-pinch  has  been  simulated  to  study 
rotational  instabilities  driven  by  centrifugal  forces.  Unstable 
m=2  modes  have  been  observed  when  the  magnitude  of  the  ion  ro¬ 
tational  frequency  exceeds  that  of  the  difference  between  the 
ion  and  electron  rotational  frequencies,  as  expected.  Pre¬ 
liminary  results  indicate  linear  growth  rates  comparable  to 
those  predicted  by  finite  Larmor  radius  MHD  and  Vlasov-fluid 
models. 


Work  supported  by  ONR  Contract  No.  N000l4-17-C05?8 
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ORBIT-AVERAGED  AND  IMPLICIT  PARTICLE  CODES 

R.  P.  Freis,  B.  I.  Cohen,  and  J.  A.  Byers 
Lawrence  Livermore  National  Laboratory,  University  of  California 
Livermore,  California  94550 
and  V.  A.  Thomas,  University  of  California  at  Berkeley 
Berkeley,  California  94720 

ABSTRACT 

Significant  progress  has  recently  been  made  in  extending  the 
applicability  of  particle  codes  to  long  time  scale  plasma  phenomena  .  This 
has  been  accomplished  by  means  of  orbit-averaging  and  implicit  techniques  that 
relax  the  usual  time  step  constraints  required  for  numerical  stability  of 
conventional  explicit  differencing  schemes.  The  important  features  of  the 
orbit-averaged  implicit  methods  are  substantial  reduction  of  particles, 
unconditional  stability  of  the  field  equation(s)  solution  at  large  time  step, 
and  tailoring  and  splitting  of  particle  time  scales  by  species  and  even  by 
velocity  class  to  accurately  and  efficiently  resolve  particle  trajectories. 
Some  of  the  more  important  applications  we  anticipate  are  two-dimensional  and 
limited  three-dimensional  kinetic  simulations  of  equilibrium  formation  and 
transport  in  tandem  mirrors  with  thermal  barriers  and  gun-produced 
field-reversed  mirrors  and  of  low  frequency  drift  waves  or  MHD  instabilities 
with  finite  Larmor  radius  effects  naturally  included.  The  new  simulation 
techniques  allow  use  of  much  more  realistic  simulation  parameters. 

We  report  the  successful  implementation  of  new  implicit  methods  and 
merging  of  orbit  averaging*  with  the  moment  equation  field  implicit  method 
independently  developed  by  Mason2  and  Denavit^  and  with  the  direct 
implicit  particle  method  due  to  Langdon,  Friedman  and  Cohen.  We  have 
analytically  determined  rudimentary  stability  and  filtering  characteristics 
and  confirmed  stable  code  operation  with  a  linearized  electrostatic  code  and 
nonlinear  one-dimensional  electrostatic  and  two-dimensional  magneto-inductive 
codes.  Orbit  averaging  by  itself  may  lead  to  numerical  instability  depending 
on  the  underlying  physics  model.  Merging  orbit  averaging  with  one  of  the 
implicit  methods  may  not  give  a  stable  or  convergent  scheme  unless  sufficient 
implicitness  is  preserved.  We  will  exhibit  code  results  on  the  simulation  of 
ion  acoustic  waves  and  the  multiple  time-scale  simulation  of  mirror  machine 
build-up,  equilibrium  and  transport. 

We  have  also  developed  a  particular  form  of  the  Mason  implicit 
technique^  in  a  linearized  particle  code  that  accurately  reproduces  growth 
rates  and  real  frequencies  for  ion  cyclotron  instabilities  and  low  frequency 
gravity-driven  interchange  modes  with  FLR  stabilization  effects  properly 
calculated.  This  code  allows  any  value  of  Wp^At  »  1. 

1.  B.  I.  Cohen,  T.  A.  Brengle,  D.  B.  Conley,  and  R.  P.  Freis,  J.  Comp. 

Phys.  38,  45  (1980). 

2.  R.  J.  Hason,  "Moment  Equation  Field  Implicit  Particle  Simulation  of 
Plasmas,"  Los  Alamos  Scientific  Laboratory  Report  LA-UR-80-2171, 

August  1,  1980. 

3.  J.  Denavit,  "Time  Filtering  Particle  Simulations  with  uugAt  »  1," 
Lawrence  Livermore  National  Laboratory  Report  UCRL-85097,  October  1980. 

♦This  work  was  supported  by  the  U.S.  Department  of  Energy  under  Contract  No. 
W-7405-ENG-48  at  Lawrence  Livermore  National  Laboratory  and  DOE  Contract 
AS03-76SF00034 -0 E- AT03 - 76 ET53064  at  U.C.  Berkeley. 
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